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ABSTRACT 

We have carried out two-dimensional, axisymmetric, hydrodynamic numerical modelling of 
the evolution of radio galaxy lobes. The emphasis of our work is on including realistic hot- 
gas environments in the simulations and on establishing what properties of the resulting radio 
lobes are independent of the choice of environmental properties and of other features of the 
models such as the initial jet Mach number. The simulated jet power we use is chosen so that 
we expect the inner parts of the lobes to come into pressure balance with the external medium 
on large scales; we show that this leads to the expected departure from self-similarity and the 
formation of characteristic central structures in the hot external medium. The work done by 
the expanding radio lobes on the external hot gas is roughly equal to the energy stored in the 
lobes for all our simulations once the lobes are well established. We show that the external 
pressure at the lobe midpoint is a reasonable estimate of the internal (lobe) pressure, with only 
a weak dependence on the environmental parameters: on the other hand, the predicted radio 
emission from a source of a given physical size has a comparatively strong dependence on 
the environment in which the lobe resides, introducing an order of magnitude of scatter into 
the jet power versus radio luminosity relationship. X-ray surface brightness and temperature 
visualizations of our simulations bear a striking resemblance to observations of some well- 
studied radio galaxies. 
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1 INTRODUCTION 

Ever since the establishment of the standard ('beam') model for 
powerful double radio sources ( |Scheuer||1974||Brandford & Rees| 
11974) it has been clear that their dynamics must depend strongly on 
their environment. Scheuer (1974) proposed two limiting cases. In 
his model A, the transverse expansion of the radio galaxy ignores 
the external medium; this would be true in practise if radio galaxies 
were strongly overpressured with respect to their external medium 
at all times in their lifetime, in which case the lobe pressure would 
drive a supersonic transverse expansion and we would expect to 
see an elliptical bow shock surrounding the sources at all times. 
In model C, on the other hand, the pressure in the lobes becomes 
comparable to the external pressure, and then, as Scheuer put it, 
"the outer parts of the cavity swell at the expense of the parts near- 
est the massive nucleus, where the thermal gas pressure is higher". 
Such a model, as Scheuer realised, is significantly harder to deal 
with analytically. 

For this reason popular analytic models of radio lobe dynam- 
ics have tended to be descendants of Scheuer' s model A. Begel-| 
[man & Cioffi (198 9) showed that it was possible to construct self- 
consistent models in which the lobes remained overpressured, and 



therefore supersonically expanding, throughout a significant part of 
their lifetimes, while allowing the jet to be confined by the mate- 
rial inside the cocoon. Models of this sort also provide the basis 
for the important and widely used work of |Kaiser & Alexander] 
(1997) (hereafter KA) who constructed an analytic model for the 
growth of a radio source in a power-law atmosphere using the as- 
sumption of self-similar expansion, which restricts the applicabil- 
ity of the model of the model to the strongly overpressured phase. 
More recently, [Hardcastle & WorraTl ( 2000) showed that the min- 
imum pressures in radio lobes were typically comparable to, or 
less than, the external thermal pressures in the centres of the host 
groups or clusters, and argued that the KA model, applied to ob- 
served extended 3CRR classical double (Fanaroff & Riley 1 974| 
class II, hereafter FRII) sources in group or cluster environments, 
required much lower cluster temperatures for self-consistency than 
are actually plausible. Since then, the routine use of X-ray inverse- 
Compton pressure measurements combined with the information 
provided by X-ray observations on the external medium (e.g. |Hard^ 
[castle et al.l|2002a|[Croston et al.||2004| ) have shown that our best 
estimates of the lobe pressures are comparable to the external pres- 
sure on scales comparable to those of the midpoints of the lobes. 
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and therefore in general significantly less than the pressure at the 
centre of the cluster or group environment; either these best es- 
timates of pressure are wrong (for example, because the contri- 
butions of non-radiating particles to the lobe pressure has been 
severely underestimated), or lobes are not strongly overpressured 
with respect to their environments, and models such as those of 
KA cannot consistently describe them on the largest scales. Ana- 
lytical calculations of the structure of FRII radio sources approach- 
ing pressure equilibrium have been presented by Alexa nder |p002] ); 
these show that in atmospheres in which pressure/density decreases 
steeply with radius, a dense layer of ambient gas is predicted next to 
the radio lobes. Its gravity squeezes the lobes and causes the onset 
of non- self- similar evolution. 

While analytical work is enlightening regarding the general 
structure and basic physics of those sources, it relies on many as- 
sumptions, e.g. about the geometry of the sources. Numerical mod- 
elling, in principle, frees us of the limiting assumptions of analytic 
models and allows us to model radio galaxies whose behaviour is 
in better agreement with observation. However, the high-resolution 
simulations required, with their large demands on simulation time, 
have typically required some compromises. Early numerical mod- 
elling (e.g. |Norman etaLl[T982l|Kossl & Mmier][T988lpndetaLl 
|1989| ) assumed a uniform-density environment, which is quite pos- 
sibly not a bad assumption for the early phases of a radio galaxy's 
growth but certainly not valid on scales of tens to hundreds of kpc. 
This was followed by 2D and 3D hydrodynamical modelling of 
sources in semi-realistic (/3-model) environments (e.g. Reynolds' 
etal. 2002 , Basson "&ATexa nder 2003 , Zan ni et al.| [2003 , Krause 
2005|. These generally showed deviations from self- similarity in 
the sense that the axial ratios of the simulated sources did not re- 
main constant, but instead grew with time. More recently, [Gaibler , | 
[Krause, & Camenzind| p009| have demonstrated in a magnetized 
jet simulation that self- similarity breaks down when a radio source 
gets close to pressure equilibrium. Qualitatively, similar results are 
seen in other MHD simulations with realistic environments (e.g. 
lO'Neill et al.] [20051 [O^ill & Jones|[20T0l ); environments derived 
from simulations of dynamically active clusters are used (e.g. |Heinz| 
|et al.|[2006 ) and the culmination of this type of work is represented 
by, for example, [Mendygral, Jones, & Dolag| ( |2012| ) whose simula- 
tions are fully 3D, contain magnetic fields, take account of electron 
transport and loss processes, and embed the radio galaxy in an en- 
vironment that is itself extracted from a cosmological simulation. 
Because of the limitations in CPU time, though, the vast major- 
ity of the work carried out in this area has involved small numbers 
of simulations considering one or at most a few different environ- 
ments. It is thus difficult to get an overview of the effects that the 
known range of environments, even simply considering such basic 
properties as scale lengths and large-radius power-law indices, has 
on the dynamics of the radio sources. 

Comparison between analytical and numerical work has also 
given rise to some puzzling results. As noted above, KA presented 
convincing physical arguments for the existence of a self-similar 
expansion law in the overpressured phase. Yet this self-similar ex- 
pansion is not easily seen in simulations (e.g. [Carvalho & O'Dea] 
[2002 V This difficulty is related to the fact that simulations often 
use a jet that is already collimated as a boundary condition to the 
simulation. But, as KA showed, the self-similar expansion phase is 
linked in a crucial way to the self-collimation of the jet by the over- 
pressured lobes (compare also model B of |Scheuer|1974| ). Thus, if a 
jet is injected into a computational domain in a parallel, collimated 
way, the jet radius cannot adapt in the proper way, which results 
in an expansion law that is in detail unphysical. For FRII jet simu- 



lations, 'KomissaroT'&Faile^ have shown that a self-similar 

scaling is indeed achieved if the jet is injected conically. They were 
even able to account in detail for the shock structure observed in 
the Cygnus A jet, produced by the self-collimation, albeit with a 
worryingly large initial opening angle. [Krause et al.| ( [2012| ) have 
re-analysed the problem using similar simulations, and found that 
the problem was related to a matter of definitions: when the same 
definitions are applied to simulations and observations, they are en- 
tirely consistent. In the present paper, we aim to set up FRII radio 
sources where the jet is self-collimated, so that self- similar expan- 
sion is in principle possible, and then to investigate their behaviour 
when they come into pressure equilibrium. 

The desire to know exactly how radio galaxies evolve in re- 
alistic environments is not a purely abstract one. Radio galaxies in 
general provide a unique way to couple the high power of an active 
nucleus (AGN) to its environment on scales of hundreds of kpc to 
Mpc, where the effects of AGN radiation are expected to be neg- 
ligible. They therefore play an important role in models of galaxy 
formation and evolution, where they provide so-called 'feedback', 
helping to prevent the cooling of large-scale gas and the consequent 
growth of the host galaxies (e.g. [Croton et al.[ [200^. However, 
whether the physics of radio galaxies is consistent with the role they 
are required to play in these models is still an open question (see 
e.g. [McNamara & Nulsen|[20T2l for a recent review). For example, 
it is not remotely clear whether the low-power, Fanaroff-Riley class 
I (FRI) radio galaxies with extended twin jets, seen in many rich 
group and poor cluster environments, are having any significant en- 
ergetic impact on the gas with short cooling times, which exists on 
scales of a few kpc from the nucleus (e.g. Har dcastle et al.[[2002b[ 
[Jetha et al.|[2007| l. Such systems are very difficult to model in detail 
numerically (though see, e.g., the work of Perucho & Martf|2007| ). 
Classical double radio galaxies, although rarer, are important in this 
context simply because they are required to have a strong effect on 
the external medium through the driving of strong shocks, which 
increases the entropy of the external medium and offsets cooling, 
and does so through a comparatively large volume of the environ- 
ment (much larger than the observed radio lobes); they are also 
simpler to model. Indeed, some of the early work on numerical 
modelling with realistic environments (e. g. [Basson & Alexander[ 
^003, Krause 2005, Zanmetal. 2005, O'Nei ll et al.[[2005| ) was 
done with the explicit aim of addressing questions about the en- 
ergetic effects of a powerful radio source. In this context, [Zanni[ 
|et al.| ( |2005 ) have shown that energy input by jets may reconcile 
cluster atmospheres from cosmological simulations with observed 
X-ray scaling relations for groups and clusters of galaxies. How- 
ever, again, the small numbers of simulations carried out in these 
studies mean that it is hard to get an overview of the differences in 
energy transfer efficiency, if any, that might be expected for differ- 
ent environments. Consequently observers, rather than making use 
of the results of numerical modelling, have a tendency to use naive 
estimates of the work done on the external medium (e.g. pdV for 
some pressure p) together with estimates of the jet power based on 
analytical rather than numerical work. 

In the present paper we try to address these disconnections 
between numerical modelling and observation. We use the very 
large amount of computational power that has become available 
over the course of the past decade, not to include still more physics 
in a simulation of an individual radio galaxy, but to carry out sim- 
ple (though high-resolution and large-volume) simulations of radio 
galaxies in a wide range of environments. Our objective in doing 
so is to draw conclusions about the dynamics and energetic input 



© 0000 RAS, MNRAS 000, 000-000 



Numerical modelling of radio galaxy lobes 3 



of radio galaxies that may be generally applicable, and to test the 
validity of some assumptions commonly used by observers. 



2 THE SIMULATION SETUP 

Our modelling in this paper uses the freely available code PLUTCQ 
version 3.1.1, described by Mignone et al.| (|2007|). We chose to 
use PLUTO because of its good handling of high Mach number 
flows and the simplicity with which new problems can be defined. 
In particular, we employed Pluto's 'two-shock', exact Riemann 
solver, third order interpolation with the van Leer flux limiter and 
third-order Runge-Kutta time- stepping. The global accuracy of the 
method is, however, second order due to the way the fluxes are cal- 
culated. 

Our general approach follows that of |Krause et al.| ( [20T2] ); that 
is, we inject two, oppositely directed, initially conical flows with 
high Mach number at the centre of the simulation, and allow them 
to recollimate in the external pressure. In a uniform-density en- 
vironment, we would expect this recollimation to take place on a 
physical scale comparable to the inner scale defined by | Alexander] 
( [2506] ), 



1/2 



where Qo is the jet power, px the external density and vj the jet 
speed. In such an environment, we would also expect the lobes to 
come into pressure balance on scales comparable to (Komissarov] 
|&Falle,.1998j ) 



Qo 

PxCl 



1/2 



where Cx is the external sound speed. As L2/L1 — 
where M is the Mach number of the jet, it is immediately clear 
that high-resolution simulations are necessary to model high-Mach- 
number (and therefore realistic) jets out to the scale on which the 
lobes are likely to come close to pressure equilibrium; the cell size 
of the simulation needs to be ^ Li and the outer radius of the 
simulation needs to be of order a few times L2. For this reason, 
we chose to reduce the complexity of the problem by carrying out 
2D (except as described below), axisymmetric, hydrodynamic sim- 
ulations. The co-ordinates are spherical polars and the jets are im- 
plemented as a boundary condition at the inner radius, injecting a 
conical flow with opening angle of 15° [chosen to ensure FRII-like 
structure: see |Krause et al.|p012| ) for a discussion of the effect of 
opening angle on the lobe morphology]; that is, at r = rinner, for 
all points with < 15° or ^ > 345°, we set p = po, Vr = Mcs, 
ve = 0, while for all other values of we use a reflection boundary 
condition at rinner- A conserved tracer quantity is injected with the 
jets and we use non-zero values of this quantity to define the lobes 
in the simulations, as described in more detail below. 

An important feature of our modelling is that the environments 
are not uniform. Instead, we represent the environment by a stan- 
dard isothermal /3 model (King profile). 



n = no 



1 + 



-3^/2 



^ http://plutocode.ph.unito.it/ 
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This introduces a new scale, rc (the core radius) into the prob- 
lem, meaning that self-similar analytic approaches (e.g. KA) are 
no longer applicable, but it allows us to represent a range of realis- 
tic environments from groups to clusters of galaxies. We impose a 
gravitational potential, as described by |Krause| ( [2005| , to stabilize 
the /3-model atmosphere. 

Within the simulations, we take the unit of length to be Li , the 
unit of density to be po = noprup, where p is the mean mass per 
particle, and the unit of speed to be Cs, the speed of sound in the 
undisturbed ambient medium. In these units, the central pressure 
Po is 1/7, or 3/5. In order to choose sensible values for rc in these 
units, we need to decide on an external calibration for these quan- 
tities. Since our interest is in low-power FRII radio galaxies in rich 
group/poor cluster environments, we adopt Q — 10^^ W (10^^ erg 



"^), a plausible value for a low-power FRII jet (e.g. Daly et al 



|2012| ). For the basic environmental properties we use kT = 2 keV, 
and no = 3 X 10^ ni~^, which are reasonable values for the type 
of environment we aim to study; we take p — 0.6. With these val- 
ues we have Cs — 730 km s~^, and L2 ^ 100 kpc; as the core 
radii of groups or clusters are of the order of tens to hundreds of 
kpc, this means that the scales that we will be studying are not too 
disparate. The choice of a reasonable Mach number is then set by 
the requirement that L2/L1 should not be too large for simulations 
at reasonable resolution. This requires a compromise regarding the 
Mach number. In reality, the Mach number should be (at least) sev- 
eral hundred, which leads to a scale separation of L2/L1 > 1000. 
At this point, we are only able to simulate jets out to several hun- 
dred Li, requiring < 50 to allow us to reach scales of a few 
times L2. In what follows we adopt a range of Mach numbers, but 
our fiducial value isM = 25, which for the parameters above gives 
Li = 2.1 kpc. Our typical simulation has a resolution of 2000 (r) 
X 1600 (0) and an outer radius of 150Li, which for the parameters 
above is ~ 3L2 ~ 300 kpc; thus we can both sample on scales 
significantly less than Li (which is necessary for the simulations 
to produce a self-collimating jet) and go out to scales significantly 
larger than L2. For these parameters, the simulation time unit is 
T = Li/cs = 2.9 X 10^ years, and we expect to simulate the 
growth of the radio source out to the edge of the computational 
volume, which (since we expect the radio source growth to be and 
remain supersonic) implies simulations lasting at least several tens 
of simulation time units, or of the order 10^ years; this is compa- 
rable to plausible maximum ages for real radio galaxies [e.g. to the 
maximum spectral ages inferred for giant radio galaxies by Jam-| 
|rozy et alT] ( |2008| ); note that spectral ages tend to underestimate the 
true dynamical age ( |Eilek|p^96t ]. 

Our choices of (3 and rc for the modelled atmosphere should 
then reflect the variety of environments observed in real groups 
or clusters. Observationally, /3 varies between about 0.3 and 1.0 
in groups hosting radio galaxies (see, e.g., |Croston et al.||2008| ); 
note that we are only simulating the large-scale environment and 
so comparisons should be made with /3 rather than /3in in this pa- 
per - the small-scale dense cores seen in observations would only 
affect the dynamics of the radio galaxy over a comparatively short 
period). Core radii span the range of tens to hundreds of kpc, as 
noted above. As described in more detail in the next section, we 
carry out simulations using a range of /3 and core radius values 
to allow us to investigate which properties of the simulated radio 
galaxy are dependent on, and which are independent of, the choice 
of environment. 

We choose to model bipolar jets because this avoids any un- 
certainty in setting the correct boundary conditions on the plane 
= 7v/2 radians (the plane passing through the origin and perpen- 
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dicular to the jet axis). To break the symmetry, we impose a sUght 
sinusoidal variation on the Mach number of the two jets, such that 
M = Mq[1 + Ashi{ujt + 0)], where t is the simulation time in 
internal units, cj and A are constants {A <^ 1), but (j) is different for 
the two jets. As the rate of energy input into the simulations goes 
as A^^, this in principle affects the calculation of the appropriate 
length scales given above, but the time-averaged energy input scal- 
ing factor is 1 + so for small A the effect is negligible. In 
all our simulations we used u — 2, A — 0.1; = for the jet 
propagating in the ^ = direction and 1 for the jet in the ^ = tt 
direction. This process has the advantage that we have two quasi- 
independent jet simulations for each run we carry out. The results 
of these may be presented separately to assess the scatter in any 
measured quantity, or may be averaged to reduce it. 

Important aspects of the astrophysics of the interactions of 
radio galaxies are omitted by this simulation setup. For example, 
there is no radiative cooling, either in the external gas or in the 
material inside the lobes. There is no particle acceleration; the ma- 
terial injected in the jets is, and remains, ordinary gas with a non- 
relativistic equation of state. There are no bulk relativistic motions 
(for example, M — 2^ corresponds to 0.06c for the parameters de- 
scribed above). Magnetic fields are not simulated, and, of course, 
the simulations are axisymmetric and two-dimensional. Our hope 
is, however, that the global dynamical properties of the sources and 
their energetic impact should be largely unaffected by these simpli- 
fications, which mainly affect the detailed behaviour of the plasma 
inside the lobes. We will discuss the effects that these limitations 
have on the interpretation of our results, together with inferences 
that can be drawn from other authors' more sophisticated simula- 
tions, in later sections. 



3 THE SIMULATION RUNS 

We carried out several basic types of simulation run: 

(i) Runs with = 25 and outer radius 150Li as described in 
the previous section. These were our standard runs and sampled a 
4x3 grid in /3 and rc, where /3 could take the values 0.35, 0.55, 
0.75 or 0.90 and rc the values 20, 30 or 40Li, corresponding to 42, 
63 or 84 kpc respectively. 

(ii) Larger runs with = 25 and outer radius 250Li. These 
were intended to allow us to investigate the late-time behaviour of 
the sources, and correspond to total source sizes of 500Li or 1 
Mpc, i.e. to scales comparable to the largest known radio galaxies. 
Because they required longer run times, we did not run these over 
the full grid of (3 and rc values. 

(iii) Runs with M / 25. These enable us to search for any 
trends with Mach number. To allow like-for-like comparison be- 
tween these and the other runs, we scaled the resolution, the core 
radius and the outer radius so as to retain (1) the same resolution as 
a fraction of Li, (2) a core radius in physical units that is matched 
to some of the = 25 simulations, and (3) the same physical 
outer radius. We chose to run all of these simulations for a 'repre- 
sentative' set of cluster parameters, /3 = 0.75 and rc = 63 kpc, i.e. 
30Li for M. — 25. For ease of comparison between simulations of 
different Mach number, all simulation lengths in plots are scaled to 
physical units (kpc) by multiplying by Li in what follows. 

(iv) 3D runs. These were designed to check that our results are 
not simply a result of the use of 2D simulations by having a num- 
ber > 1 of cells in the <j) axis of spherical polars. This makes the 
problem computationally considerably harder, since we need to re- 
tain the resolution of the 2D runs in the r and directions in or- 



der to sample the jet injection region appropriately, as discussed 
above. Consequently, the final 3D runs we present have compara- 
tively coarse resolution in the axis, and are also restricted to one- 
sided jets. They are roughly matched to two of the 2D runs (one has 
a smaller outer radius) and have A^ = 25. The axisymmetry of the 
problem is broken in a similar way to the approach taken with two- 
sided jets, by adding a small precessing sinusoidal perturbation of 
the injection Mach number as a function of time and (j). We only 
use these runs in comparisons with the corresponding 2D runs. 

Table [T] summarizes the runs that we refer to in the rest of 
the paper. Each run is given a short code (e.g. 'M25-90-30', corre- 
sponding to A^ = 25, /3 = 0.9 and rc = 30) which indicates the 
key parameters of the run, and these codes will be used, for brevity, 
throughout the remainder of the paper. 

All our production runs were carried out using the STRI clus- 
ter of the University of Hertfordshire, using either 128 or 192 Xeon- 
based cores, typically for about 24-48 hours per run. Some early 
runs to explore parameter space were carried out on smaller num- 
bers of cores either on the STRI cluster or on the UCL supercom- 
puter Legion (access provided by the Miracle project). 

Post-processing was also carried out on the STRI cluster. 
PLUTO was configured to write out the complete state of the sim- 
ulation every 0.1 simulation time units, and these images of the 
simulation grid were used to compute derived quantities such as 
the dimensions of the lobe, the energy stored in the lobe and in 
the shocked region, and so forth. These quantities are used in the 
following section. One crucial choice concerns the exact definition 
of the lobe: we defined material as being inside the lobe if it was 
within the outermost surface where the lobe tracer quantity was 
> 10 We use a value < 1 here to account for mixing, which 
leads to dilution of the original jet material. The exact results are 
insensitive to the value used as long as it is not too close to unity, 
which would lead to spurious 'heating' of the external material due 
to mixing with lobe material, a process which does not appear to 
occur in real radio galaxies, where a sharp boundary between the 
lobes and the external medium is always seen. Material is taken to 
be within the shocked region if it is not in the lobe and if it is inside 
the outermost surface where the radial velocity > 10~^ simulation 
units. Calculation of energies for the shocked region take as their 
zero point at any given time t the energy stored within the boundary 
of the shock at time t in the image from t — 0, so that we ignore 
the pre-existing internal energy of the environment. 



4 RESULTS 

4.1 General properties of the simulations 

All simulations produce reasonably plausible 'radio lobes' after the 
initial establishment and collimation of a jet. Figs [T] and [2] show 
maps of the density and temperature as a function of position in 
Cartesian co-ordinates for several different times in an example 
simulation (M25-35-40), illustrating the general trends we observe. 

At the start of the simulations the jets initially form quite long, 
thin lobes, but these then expand transversely as the material that 
has travelled up the jet is thermalized by shocks driven into the 
lobes by new jet material. This initial phase is a necessary conse- 
quence of an initially overdense conically expanding jet: the width 
of the lobes is related to the ratio of the jet density to the ambi- 
ent density (e.g.'Kra usel|2003| », and only very underdense jets form 
wide radio lobes. The thin lobe phase that we see is therefore a re- 
alistic feature, but it appears on unrealistically large scales, set by 
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Table 1. The parameters of the simulation runs discussed in the paper. The columns give the type of simulation, as discussed in the text (Section [Sj; the 
sidedness of the jet (1 for 3D runs, 2 for all others); the Mach number and the corresponding simulation length and timescale, as discussed in the text (Section 
|2j; the outer radius of the simulation, in internal units; the resolution; the properties of the environment; the timescale for the simulation to run to completion, 
in internal units; and the code used to designate the simulation in the rest of the paper. Where a column is blank, it has the same value as the nearest non-blank 
entry immediately above it. 



Run type Sided 


M 


Li 


r 


Outer radius 


Resolution 


/5 


rc 


^end 


Code 






(kpc) 


(Myr) 




(cells, r X 0{x(j))) 






(r) 




Standard 2 


25 


2.1 


2.9 


150 


2000 X 1600 


0.35 


20 


44.2 


M25-35-20 














0.35 


30 


51.7 


M25-35-30 














0.35 


40 


60.0 


M25-35-40 














0.55 


20 


31.3 


M25-55-20 














0.55 


30 


29.2 


M25-55-30 














0.55 


40 


48.1 


M25-55-40 














0.75 


20 


22.9 


M25-75-20 














0.75 


30 


37.7 


M25-75-30 














0.75 


40 


49.6 


M25-75-40 














0.90 


20 


21.7 


JV125 -90-20 














0.90 


30 


34.1 


M25-90-30 














0.90 


40 


42.6 


M25-90-40 


Large 2 


25 


2.1 


2.9 


250 


3000 X 2400 


0.35 


40 


95.0 


M25-35-40-L 














0.55 


40 


77.1 


M25-55-40-L 














0.75 


40 


60.1 


M25-75-40-L 














0.90 


40 


45.0 


M25-90-40-L 


Varying M 2 


10 


8.5 


11.4 


38 


500 X 400 


0.75 


7.59 


7.6 


MlO-75-30 




15 


4.6 


6.2 


70 


930 X 740 


0.75 


13.94 


18.2 


M15-75-30 




20 


3.0 


4.0 


107 


1430 X 1150 


0.75 


21.47 


18.5 


M20-75-30 




30 


1.6 


2.2 


197 


2630 X 2100 


0.75 


39.44 


46.7 


M30-75-30 




35 


1.3 


1.7 


248 


3310 X 2650 


0.75 


49.70 


47.3 


M35-75-30 




40 


1.1 
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the choice of our Mach number (compare above); on these scales, 
it is an artefact of our particular approach. This feature appears re- 
peatedly in the analysis below. 

Very high temperatures are observed within the lobes at all 
times, as expected given the high Mach number of the jet termi- 
nation shock. However, a stable terminal shock (corresponding to 
the hotspot of FRII radio galaxies) does not develop; the jet tends 
to terminate soon after entering the lobes (this can be seen in the 
temperature maps - the unshocked jet is cold due to adiabatic ex- 
pansion) and is unstable, with repeated jet terminations at different 
positions inside the lobe driving shocks into the lobe material itself. 
We do not regard this as a major problem, because the object of our 
simulations is to consider the dynamics of an overpressured lobe, 
rather than to reproduce the detailed observational properties of 
FRII radio galaxies. The density plots also show quite clearly that, 
on scales comparable to L2, the lobes start to be driven away from 
the central parts of the (shocked) atmosphere. By the late stages 
of the simulation, there is a substantial gap between the two lobes 
and the jets are propagating through the shocked external medium 
without a cocoon to protect them. Finally, we see that the edges of 
the lobes are Kelvin-Helmholtz unstable, as is frequently observed 
in purely hydrodynamical simulations of radio lobes: the effects of 
entrainment of shocked external material (denser, colder) can be 
seen in the density and temperature plots. We do not believe that 
this represents the behaviour of real radio sources, noting that mag- 
netic fields would suppress the K-H instability ( |Gaibler et al.][2009l ) 
but again we do not expect the K-H instabilities to have a signifi- 
cant effect on the lobe dynamics, which are governed by large-scale 
pressure and momentum flux balance between the lobe and the ex- 
ternal medium, though the instabilities do clearly have an effect on 
the detailed appearance of the simulated sources. Only if the K-H 



instabilities grew to size scales comparable to those of the lobes 
would we expect the overall dynamics to be affected, and this does 
not occur in our simulations. 

Turning to the effect on the external medium, we see that ini- 
tially the lobes drive the expected shell of hot, shocked gas, though 
even at early times high temperatures are only seen close to the 
ends of the radio lobes. The transverse expansion of the lobe is 
only trans-sonic from early times and so the post-shock gas away 
from the ends of the lobes is free to expand almost adiabatically, 
driving at most a weak shock into the external medium (in the ex- 
ample shown, we see that the shocked gas extends to a radius of 
60 simulation units in the transverse direction at t = 50). At late 
times, although there are temperature and density fluctuations in 
the shocked medium (due to weak shocks related to vortices in the 
lobes, which are probably not realistic in detail), its overall temper- 
ature and density away from the tips of the lobes is not markedly 
different from the initial conditions, except of course where the 
lobes have formed cavities. Note in particular the formation of a 
bar-like structure and arms of colder material around the base of 
the lobes by t = 50, structures whose visibility would presumably 
only be enhanced if radiative cooling were included in our simula- 
tions. We will return to the observational consequences of this in 
Sectionl431 

4.2 Dynamics of the radio lobes 

In Fig s [3] and |4] we plot the mean length and volume of the lobes in 
the M — 2^ runs as a function of simulation time. Here, for each 
simulation run, we have taken the average length and volume of the 
two lobes in order to reduce the scatter and highlight the interest- 
ing features. Considering the linear growth of the lobes first (Fig. [3]) 
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Figure 1. Density in the simulation M25-35-40 for simulation times t = 10, 20, 30, 40 and 50. The images show a slice through the midplane of the notional 
3D volume; as the simulations are axisymmetric, all images are reflection- symmetric about the jet axis. Colours show a logarithmic scale of density between 
0.025 and 6.3 simulation units, with dark blue being densest and red least dense. Black and white contours show the lobe and shocked region boundaries 
respectively, as defined in Section[3] Units of the axes are simulation units (Li). 
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Figure 2. Temperature in the simulation M25-35-40 for simulation times t = 10, 20, 30, 40 and 50. Notes as for previous figure, but colours show temperatures 
between 0.32 and 3.2 simulation units, with dark blue being hottest and red coldest. 
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Figure 3. Mean lobe length as a function of simulation time in = 25 simulations for (top left) ^ = 0.35, (top right) /3 = 0.55, (bottom left) /3 = 0.75 
and (bottom right) (3 = 0.90. Also plotted ('KA'), with arbitrary normahzation, are the |Kaiser & AlexanderH1997| predictions for a power-law atmosphere. 



we can generally pick out three distinct slopes in the plots: an initial 
rapid growth (where the jet grows quickly without forming much in 
the way of lobes, as we noted above), followed by a phase of slower 
growth, followed by more rapid expansion again. As the transition 
between these last two phases happens on scales comparable to rc 
(and in fact clearly happens earlier for the sources with smaller 
Tc; see Fig. [3]), it is natural to assume that it is the result of the 
lobes emerging from the approximately uniform-density material 
within rc and starting to probe the approximately power-law den- 
sity regime seen on large scales. Consistent with this, we have plot- 
ted the functional form of the length against time, Ll oc t^/(^-^^\ 
from the analytic modelling of KA on these plots (with arbitrary 
normalization) and we see that this is reasonably comparable to the 



slope at late times - this is particularly clear in the four runs that 
go out to radii of 250 simulation units. (We note that these four 
runs seem systematically slightly offset from the equivalent runs 
with outer radii of 150 units - we attribute this to the slightly lower 
resolution of the larger runs, which seems to mean that the lobe is 
slightly longer at all simulation times shown on these plots.) Thus 
we reproduce not just the qualitatively expected behaviour - higher 
/3 should mean faster growth at late times - but also do surprisingly 
well at reproducing the quantitative predictions of KA. 

In the volume plots (Fig.[4| we see a very similar three-phase 
picture: the lobe expands with relatively slow growth of volume 
at early times (we explore the meaning of this below), then grows 
rapidly after about t = 2 in most simulations, and finally expands 
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Figure 4. Mean lobe volume as a function of simulation time in = 25 simulations for (top left) /3 = 0.35, (top right) /3 = 0.55, (bottom left) /3 = 0.75 
and (bottom right) (3 = 0.90. Also plotted ('KA'), with arbitrary normalization, are the jKaiser & AlexanderH1997^ predictions for a power-law atmosphere. 



more slowly and with roughly a power-law dependence on t. The 
power law we see, though, is systematically flatter than the predic- 
tion from KA, Vl oc L| oc t^/^^-^^\ KA's prediction of course 
comes directly from their assumption of self- similarity in the lobes, 
so the fact that Vl does not appear to scale as L| points to the 
fact, already seen in Section [43] that these lobes do not grow self- 
similarly. 

The two 3D simulations show very similar lobe dynamics to 
their corresponding 2D simulations, as shown in Fig. [5] Over most 
times the growth of length and volume in 3D is indistinguishable 
from the 2D results. The only noteworthy difference is that the vol- 
ume seems systematically larger at very early times, particularly in 
the M25-75-20-3D simulation which has higher spatial resolution. 



This reinforces the impression that the initial rapid growth of the 
lobe is not well modelled, since the results depend on the numeri- 
cal approach taken. However, the 3D and 2D simulations agree very 
well after this initial phase. 

The non- self- similar growth of the lobes is illustrated in an- 
other way by Fig. [6] (top panel) which shows the axial ratios as a 
function of lobe length for all the = 25 sources in 2D simu- 
lations. We define the axial ratio as the ratio of the lobe length to 
the full lobe width, measured half-way along the lobe length: this 
is very similar to the definition used by |Mullin, Riley, & Hardcastle] 
( 2008). With this definition, large numbers imply thin lobes, small 
numbers correspond to fat ones. Two important features are seen 
on this plot. Firstly, we note the very strong peak in the axial ratios 
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Figure 5. Length and volume as a function of simulation time for the two 3D simulations (solid lines) and the two lobes of the corresponding 2D simulations 
(dashed lines). 



at small lobe lengths, which corresponds to the phase where the 
source is growing rapidly in length at close to constant lobe vol- 
ume. This phase, which is probably not realistic, is generally over 
by Ll ~ 50Li (or ^ 100 kpc) as the initial growth of the source 
slows down and the lobe begins to inflate. More interestingly, at 
Ll > 50Li we see a steady growth of the axial ratio as a func- 
tion of lobe length; that is, throughout most of the lobes' growth, 
the lobes are getting thinner as they expand. This is because, while 
the growth in length of the source is determined by the balance be- 
tween the lobe pressure and the momentum flux of the jet with the 
density and pressure at the end of the lobes, the transverse expan- 
sion of the lobes is driven by the internal lobe pressure only and the 
external density and pressure are higher. The change in axial ratio 
with time or lobe length is not consistent with self- similarity, but it 
is consistent with what is observed (Mullin et al.||2008). Our sim- 
ulations do not span the full range of axial ratios observed in real 
sources, but they do lie in the right region of parameter space, and 
of course a good deal of scatter is introduced into observed axial 
ratios by projection effects. 

Mach number has some effect on the axial ratio, in the sense 
that simulations with very low M do not form realistic lobes at any 
stage of the run; their lobe volumes remain low and so their axial 
ratios are high. This is probably because the material in low Mach- 
number lobes is never shock-heated to high enough temperatures 
that the pressure in the lobes (which drives the transverse expan- 
sion) becomes comparable to the jet momentum flux. Changing the 
jet properties, e.g. the opening angle, or including some jet preces- 
sion would probably allow lobes to form in some circumstances, 
but in any cases such low jet Mach numbers (for the jet powers we 
are considering) are not physically very plausible. For M > 25, 
the late-time axial ratio appears more or less independent of Mach 
number, and is within the range seen in the various A4 — 25 sim- 
ulations (though we note that the early phase of fast lobe growth 
seems even more extreme for M > 25). 

Finally, we consider pressure balance in the lobes. We com- 
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Figure 8. The ratio between mean pressure in the lobes and undisturbed 
external pressure at the lobe midpoint as a function of lobe length in kpc for 
varying Mach number. 



pare the mean pressure in the lobes (since the sound speed in the 
lobes is high, we do not expect any strong pressure gradients in- 
ternal to the lobes) to the undisturbed external pressure at the mid- 
points and endpoints of the lobes. This approach is similar to what 
has been done for radio galaxies, almost exclusively FRIIs, for 
which the mean lobe pressure can be estimated from observations 
of inverse-Compton emission and the external pressure from obser- 
vations of thermal X-ray emission (e.g. [Hardcastle et aL| |2002a[ 
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Figure 6. Axial ratio as a function of lobe length in kpc for (top) all the = 25 simulations and (bottom) simulations with varying Mach numbers but 
/3 = 0.75, rc = 62 kpc. 



ICroston et al) |2004l |Konar et al.| [2009) [Shelton et al^ [20TT] ). 

Observers use the undisturbed external pressure because this can 
easily be estimated observationally, e.g. by masking out the lobes 
and fitting spherically symmetric projected /3 models, although of 
course the lobes are not expected to be in direct contact with this 
material. What we find is quite striking (Fig. [TJ: after the initial 
expansion, and at ^ L2 the lobes come into, and remain in, 
rough pressure balance (within better than a factor 2) with the ex- 



ternal pressure at the lobe midpoint, though they are consistently 
strongly overpressured with respect to the external pressure at the 
ends of the lobes (as expected since they are seen to drive strong 
shocks at all times: Section [4T| l. Thus the simulations reproduce 
a key observational result. We also see (Fig. |8]) that this result is 
independent of Mach number so long as the Mach number is high 
enough that realistic lobes form (M > 25). The 3D simulations 
reproduce these 2D results. 
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Figure 7. The ratio between mean pressure in the lobes and undisturbed external pressure at (left) the lobe midpoint and (right) the end point of the lobes as a 
function of lobe length. Only the results from the long simulations are plotted here, but other runs are similar. 



Overall, the simulations seem to do reasonably well at late 
times at reproducing the behaviour we expect from real radio galax- 
ies. At early times, they are not realistic, as noted above. Only when 
the lobe has grown to a size comparable to the core radius do we en- 
ter a regime that is comparable to the behaviour of real radio galax- 
ies. However, once the simulations enter that regime, they seem to 
stay there, independent of jet Mach number, so long as it is high 
enough, and of the dimensionality of the simulation. This means 
that it is reasonable to consider the energetic effects of the radio 
source in its environment, as long as we confine ourselves to the 
period after t ^ 5 or so when the lobes are fully established. 



4.3 Energetics of the lobe/environment interaction 

In this subsection we consider the energetics and environmental 
impact of the lobes. 

The first and most basic test is to check whether the jets are 
carrying the expected level of energy into the simulation as a func- 
tion of time. The simulation energy unit e is poLfcl, or around 
4.6 X 10^^ J for our reference simulations with M = 25. So we 
can easily compare the rate of growth of energy in the lobes and 
environment with the expected jet power Qo, or, equivalently, we 
can calculate the expected input power in simulation units (A4^/8) 
to see whether the simulations are conserving energy as expected. 
The results of such a comparison are shown for some representa- 
tive simulations in Fig.|9] There are several points to note from this 
figure. Firstly, we see that there is a short period at the start of every 
simulation plotted where the total energy in the simulations stays 
constant, rather than increasing. In detail, it seems that the jet stops 
flowing into the simulation volume at these times: because the jet 
is implemented as a boundary condition, it is possible for the jet's 
entry on to the grid to be blocked, e.g. by backflow, with no en- 
ergetic consequences. Once the jet is established, though, we see 
an approximately linear growth of energy with time, though note 
the superposed oscillation which is a result of the slight sinusoidal 



variation imposed on the Mach number (Section [2]); we plot the 
two lobes separately in this figure so that the different phases can 
be seen. Finally, we note that although the gradients in different 
simulations are in close agreement, they disagree with the expected 
slope at around the 10 per cent level, even after a non-zero start time 
is imposed to account for the fact that a persistent jet only starts at 
around t = 3. As we see slightly different slopes for our runs with 
slightly different resolutions, we believe this discrepancy to be the 
result of resolution effects connected with the implementation of 
the jet as a boundary condition. In particular, it seems that the very 
edges of the jet have a lower density than they should have even 
at the inner radius, so that less energy enters the simulation than 
our simple calculation suggests. However, such structures at the jet 
inlet should stay constant with time, and this is consistent with the 
fact that the missing power is also constant with time. If there were 
significant errors in conservation of energy in the simulations, or 
in the calculations done in post-processing, we would expect the 
curves in plots like those of Fig.|9]to be non-linear with time, given 
that the growth of the lobes and the shocked regions is very clearly 
non-linear. As this is not the case, we are confident that the simu- 
lations can be used to study the energetic impact of the expanding 
lobes after early times; as noted above, the dynamics of the lobe are 
not realistic at early times anyway. 

Our post-processing accounts separately for the thermal, ki- 
netic and potential energies in both lobes and in the corresponding 
shocked regions. Thus we can investigate how the work done by the 
jet is distributed between these different contributions to the total 
energy. An example of this type of analysis is shown in Fig. [9] We 
see that, as might be expected, the thermal energy in the lobes and 
the shocked region dominates over the other components. Kinetic 
energy is non-negligible in both, however, accounting for some- 
thing like a third of the total. At late times the potential energy of 
the shocked material, i.e. the work done in lifting shocked material 
out of the centre of the galaxy, also becomes important. There are 
strong variations in the kinetic energy of the lobe material, which 
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Figure 9. Left panel: the total energy stored in the jets, lobes and external environment as a function of simulation time for each lobe of the four = 25 
simulations with Vc = 40, together with the expected slope as discussed in the text. Right panel: the energy budget of the lobes and shocked region in the 
simulation M25-35-40 as a function of simulation time. 



partly reflect the imposed variations in Mach number of the injected 
jet, but also show the effects of the repeated internal shocks in the 
lobes. 

Perhaps the most obvious feature of this plot is that the ther- 
mal and kinetic energy stored in the lobes keeps remarkably well in 
step with that in the shocked region over most of the lifetime of the 
source. We can quantify this for the simulated sources in general 
by plotting the ratio of the total energy stored in the shocked region 
to that in the lobes (Fig.[To]). We see that at all lengths (times) after 
the regular lobe growth is established (L = 30Li, corresponding 
to t 3), the ratio of external to internal energies is close to unity, 
with no very strong dependence on the parameters of the environ- 
ment (there does seem to be a weak dependence, in the sense that 
this ratio tends to be higher for higher /3 values). This result con- 
firms earlier findings in MHD simulations by O'Neill et al. (2005 ) 
and |Gaibler et al.||2009 ), but is based on a much more comprehen- 
sive sampling of parameter space. We return to this very important 
result in Section [5] 

Like others we have considered above, this result appears 
more or less independent of Mach number within the range we have 
used (Fig.|10|): even in the low-Mach-number systems where real- 
istic lobes are not formed, the energy ratio does not differ strongly 
from unity at late times. It is also reproduced in the 3D simulations 
(Fig[TTJ. 

Finally, we can compare the work done on the external 
medium with commonly used methods for estimating it. As we 
have seen, the internal energy in the lobes tracks that in the shocked 
region reasonably well over most of the lifetime of the source (Fig. 
|9]). Therefore, if the internal energy of the lobe can be estimated, 
the work done on the external medium can be estimated too. Many 
literature estimates of the work done (e.g. |Birzan et al.|[2004| ) as- 
sume that it is ^ pV, where p is some relevant pressure and V is 
the volume of the lobe, or of the observed X-ray cavity. This must 
be of the right order of magnitude, but setting aside the problems 
in estimating V in real radio galaxies, where projection effects are 



important, the question is what value of p to adopt. If the inter- 
nal pressure of the lobe (e.g. from inverse-Compton emission) is 
known, then this is easy, but in this case of course we know the 
internal energy density of the lobe as well. If the internal pressure 
is not known, then our results above suggest that the undisturbed 
external pressure at the lobe midpoint can give a reasonably good 
estimate of it. In this case, we can use our simulations to see how 
badly in error a simple pe^tV estimate is. 

The ratio between the work done on the external medium and 
Pe^tV is plotted in Fig. [12] If the radio source were in pressure bal- 
ance with the external medium at the midpoint, and if all the energy 
in the lobe were thermal energy, then we would expect this ratio to 
have a constant value of 1.5 for roughly equal internal and external 
energy. However, the fact that some of the work done on the lobes 
is in the form of bulk k.e. and that the midpoint pressure balance 
is only approximate causes deviations from this relationship; fac- 
tors in the range 2-10 are typical. In detail, the trends in Fig. [12] 
are readily understood. At early times, all simulated radio sources 
are strongly overpressured and drive strong shocks. Consequently, 
the ratio is high; it decreases as the source comes close to pres- 
sure balance. The simulations with smaller core radii and larger (3 
show higher ratios presumably because they are further from pres- 
sure balance over most of their length; therefore, at late times, for 
our smaller core radii and larger values of (3 we generally find val- 
ues closer to 10, whereas for large core radii and small betas the 
ratio tends towards two. To summarize, although the exact numer- 
ical factor is likely to be different for real radio galaxies, with a 
relativistic adiabatic index, it is safe to say that pe^tV pressure es- 
timates are likely to underestimate the work done on the external 
medium by a factor of a few; moreover, this factor will depend on 
the age of the radio galaxy and the environment in which it resides. 
Of course, it should be noted that our results apply principally to 
FRII radio galaxies, while many of the objects showing cavities are 
lower-power FRIs; although qualitatively we might expect similar 
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Figure 10. The ratio between the energy stored in the shocked region and the lobes as a function of lobe length; top, for all = 25 simulations; bottom, for 
simulations with different Mach number and /3 = 0.75, Vc = 62 kpc. 



caveats to apply to ape^tV analysis of FRIs, our simulations do not 
allow us to make quantitative statements about such systems. 



4.4 Radio emission 

Our simulations do not contain all the physics needed (i.e., elec- 
tron acceleration, loss processes and magnetic fields) to carry out 
a proper visualization of the synchrotron emission from the lobes. 



It is, however, interesting to ask how the overall radio luminosity 
of the lobes would evolve with time, and we can do this by assum- 
ing that the pressure that we measure in the lobes is contributed by 
electrons and magnetic field. Since the simulations do not constrain 
the contributions of these separately, we need to assume equiparti- 
tion or some constant departure from equipartition (as measured 
by inverse-Compton observations, e.g. jCroston et al. 2005 ). Let the 
energy density in electrons be Ue and that in magnetic field be Ub , 
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Figure 11. As for the top panel of Fig.[To] but the two 3D simulations are shown compared to the results for individual lobes of the comparator 2D simulations. 
Solid lines show the 3D simulations and dashed lines the 2D ones. No significant differences between the 2D and 3D results are seen after L ~ 50 kpc. 
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Figure 12. The ratio between the energy stored in the shocked region and pext ^ as a function of lobe length for the standard- sized 7W = 25 simulations. 



then this assumption corresponds to 



r]Ue = Ub 



2/io 



where 77 describes the energetic departure from equipartition, 77 < 
1. 

If we then assume a power-law electron distribution, it can be 
shown using standard textbook results (e.g. Longair 201 0) that the 
volume emissivity goes as where q is the electron energy 
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Figure 13. Radio luminosity (simulation units) as a function of mean lobe length for the = 25 simulations. 




10' 



25 



10 



100 

Lobe length (kpc) 



10( 



Figure 14. Radio luminosity (physical units) as a function of mean lobe length for the 7W = 25 simulations, converted from simulation units using the factor 
£ for ?7 = 0.1 as described in the text. The line styles correspond to the same simulations as in Fig. [13] Overplotted are the 3CRR radio galaxies above the 
nominal FRI/FRII cutoff of 5 x 10^^ W Hz~^ sr~^, where we have scaled down the largest angular size by a factor to account for the fact that the 
simulations are for a single mean lobe length and that a typical 3CRR source will be projected. Note that the different curves represent radio sources with the 
same jet power, but in different environments. The radio luminosity increases in environments with a larger core radius and shallower density decline. 
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power-law index; the luminosity can be calculated by integrating 
this over the lobe volume. Taking ^ = 2.2, corresponding to a low- 
frequency spectral index a — 0.6, the emissivity scales as p^'^ . We 
begin therefore by plotting J p^'^dV over the simulation volume 
as a function of lobe length (Fig.[T3]). These plots correspond to the 
'tracks in the P-D diagram' commonly used in studies of the evo- 
lution of radio galaxies, albeit still in simulation units at this point. 
We see first of all that the 'luminosity' in these plots varies strongly 
with source length (and therefore time). Ignoring the strong varia- 
tion before the lobe is properly established at L < 50 kpc, we 
see a general trend for the luminosities first to scale up and then 
to level out. At the longest lengths, the luminosities may level out 
or decrease slightly, the most obvious decreases being seen in en- 
vironments with high p and small rc. At L ^ 200 kpc, there is 
almost an order of magnitude difference in the radio luminosities 
of sources with, by construction, identical jet powers. We also note 
the smaller- scale, random or slightly periodic variations in the ra- 
dio power; these reflect fluctuations in the pressure in the lobes and 
should not be taken too seriously, since the detailed hydrodynamics 
in the lobes is probably not correct. 

It is important to note that the analysis described here takes no 
account of radiative losses. The approximation that the electron en- 
ergy distribution is a power law is not a bad one even when radiative 
losses are effective, because the energy density in the electrons is 
dominated by low-energy electrons with long loss timescales; but 
of course the frequency v needs to be low in order to sample the 
part of the synchrotron spectrum that is still adequately represented 
by a power law. Since we wish to compare these radio luminosi- 
ties to real luminosities, we use a notional observing frequency of 
178 MHz, which should be low enough to avoid significant spec- 
tral ageing effects in all but the oldest sources, and allows direct 
comparison to observations of the 3CRR sample and to the work of 
[Kaiser, Dennett-Thorpe, & Alexa nder (1997 ). 

Given that the simulation unit of pressure is P = ^ynokT, or 
1.6 X 10~^^ Pa, and the unit of volume is Lf , it can be shown that 
the simulation unit of radio luminosity in physical units (W Hz~^ 
sr~^) is 



C = c{q)- 



eocrrie 



q + 5 

q+1 



-Li 



where e is the charge on the electron, me its mass, c the speed of 
light, eo the permittivity of free space, /xo the permeability of free 
space, / the integral of EN{E) over the energy range of the elec- 
trons, ^min to Ema^, and c{q) a dimensionless constant, the prod- 
uct of a number of gamma functions, of order 0.05. Taking g = 2.2, 
= 178 MHz, ?7 = 0.1 (given that the measured magnetic field 
strengths can be up t o a factor of a few below equipartition, see 



Croston et aLl ( |2005| ), ^min = lOm^c^, and ^max = 10^ 



(the results are insensitive to this choice of electron energy range), 
we find that £ = 6.6 x 10^^ W Hz"^ sr"^ for the = 25 sim- 
ulations. If we take 77 = 1, the conversion factor roughly doubles, 
£ = 1.4 X 10^^ W Hz-^ sr ^, so that we are not excessively 
sensitive to equipartition assumptions. 

This conversion factor allows us to generate a plot of the evo- 
lution of radio luminosity in physical units for the = 25 sources 
(results for other Mach numbers are similar but are omitted for clar- 
ity). We overplot data for the 3CRR source^ (Fig. 14 note the 
change of x-axis scale from linear to logarithmic). We see immedi- 
ately that the luminosities are in the right general area of parameter 



space for powerful radio galaxies; of course, the 3CRR FRIIs we 
have overplotted should represent a wide range of jet powers and 
environments, down to objects with Q ^ 10^^ W in very poor 
groups (e.g. |Hardcastle et aLl|2012| ), so it is not surprising that the 
simulations do not sample the full P — D diagram. Within the lim- 
itations of the simulations and the calculation we have carried out 
(i.e., a comparison is only possible for L > 50 kpc, and we neglect 
losses, which may well be significant on scales > 300 kpc) we view 
the position of our simulation tracks on this plot as good evidence 
that we are generating physically reasonable radio sources. We note 
that the typical luminosity, ^ 10^^^ W Hz~^ sr~^, is ex actly what 



we would expect from a 10 W jet using the results of Willott et al. 



Data from http://3crr.extragalactic.info/ 



( |1999| ), which, if a coincidence, is certainly a remarkable one. 

Having said this, it is striking that the most luminous sources 
in the simulations (those located in large-rc, low-/3 environments) 
approach the highest luminosities seen for 3CRR sources, which 
is perhaps surprising given that we are only trying to simulate 
intermediate-power FRIIs with Q — 10^^ W; however, these en- 
vironments are very rich, perhaps unrealistically rich for most low- 
power 3CRR sources, and so we do not regard this as a significant 
problem. We also note that any non-radiating particle content in 
the lobe (see |Hardcastle & Croston|20T0| for a discussion of the evi- 
dence for these in FRIIs) would tend to reduce the required electron 
number densities and hence the radio luminosity. 

Finally, it should be noted that the very strong dependence of 
luminosity on environmental properties on 100-kpc scales, which is 
independent of these caveats, provides both a quantitative estimate 
of the boosting effect of rich environments discussed by |Barthel &| 
|Arnaud| ([1996) and a useful warning of the dangers of attempting 
to infer jet power from radio observations alone. On these scales, 
there is at least an order of magnitude scatter in the radio luminosity 
expected for a fixed jet power. We return to this point in Section|5] 

4.5 X-ray emission 

We can assess the observable effect of the radio galaxy's interaction 
with its environment in these models by visualizing the X-ray emis- 
sion from the environment. We choose to model X-ray emission as 
seen by Chandra: Chandra'^ angular resolution of ^ 0.5 arcsec 
gives a spatial resolution of ^ 2 kpc at around z — 0.2, so we adopt 
a spatial resolution of 1 simulation unit for the = 25 simula- 
tions. We assume that lobe material does not contribute to the X-ray 
emission, and therefore exclude it using the tracer. We use XSPEC 
with a real Chandra response file to model the conversion between 
emission measure and (on-axis) count rate in the 0.5-5.0 keV en- 
ergy band as a function of temperature for an APEC model with 
0.3 solar abundance (neglecting Galactic absorption). We can then 
compute the expected X-ray image (for infinite signal-to-noise) 
by integrating along lines of sight through the simulated volume: 
we neglect emission outside the simulation (corresponding to per- 
fect background subtraction). We can also compute the emission- 
weighted (strictly, counts-weighted) mean temperature along the 
line of sight; Chandra does not quite measure this quantity, but it 
is close enough to give an idea of the expected observational ef- 
fects. We show representations of the counts image and the map of 
emission-weighted temperature for several times during the simu- 
lation M25-35-40 and for two different angles to the line of sight in 
Figs.[l5]and[T6] 

Bearing in mind that we are 'observing' with infinite signal 
to noise and that therefore many of the details seen in the simula- 
tions are not likely to be visible in reality, these visualizations are 
qualitatively very consistent with real observations. We note first 
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Figure 15. X-ray visualization of the simulation M25-35-40 for simulation times (from top to bottom) t = 10, 20, 30, 40 and 50. The left-hand images 
show Chandra counts in the energy range 0.5-5.0 keV, with a linear transfer function. The right-hand images show the corresponding emission-weighted 
temperature, with dark blue being coldest (1.8 keV) and red being hottest (> 3 keV); see the text for discussion. The jet axis is at 90° to the Hne of sight 
and the image resolution is 1 simulation unit (~ 2 kpc). Compare Figs[^and[2] which show slices through the midplane of this simulation at high resolution. 
Co-ordinates shown are simulation units. 
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Figure 16. As for Fig. [15] but the j 



of all that, though there is an approximately elliptical shock with 
a clear surface brightness contrast at the boundary, the measured 
temperature difference across the shock is not likely to be very 
great; the emission-weighted temperature of the shocked region is 
^2.5 keV compared to the unshocked 2-keV gas. This compares 
very well to the modest temperature increases seen in this region in 
Cygnus A ( [Wilson et al.||2006] ). The regions of strong shocks and 








-150 -100 -50 50 100 150 



axis is at 45° to the line of sight. 

high temperatures (> 3 keV) are confined to the very tips of the 
lobes, and even then are not always present; it is easy to imagine 
that these high temperatures will be hard to observe in real radio 
galaxies, where measurements of spectra require integration over 
large regions. One is more likely to see moderately increased tem- 
peratures over a larger region at the head of the lobes, as suggested 
in the case of Cygnus A by pelsole & Fabian | ( [2007| and in 3C 452 
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by'Shelton et al.| ( [2011| . Cavities are visible, but again quite low 
in surface brightness contrast, and so not necessarily expected in 
real observations (particularly as we have not attempted to model 
inverse-Compton emission from the lobes, which in real observa- 
tions will fill in some of the cavity emission). The most obvious 
feature of the simulated observations, which should be present in 
any real observation if these simulations are well matched to real- 
ity, is the emergence of a bright central 'bar' of X-ray emission, 
together with 'arms' extending along the inner part of the lobes. 
This material is cool, with a temperature equal to, or even in places 
slightly less than, the original temperature; it is generated from the 
pre-existing dense core of the host cluster, which is initially heated 
by the radio-lobe shock but then cools adiabatically as the shocked 
region expands ahead of it, in the manner discussed by "Alexan-' 
|derj ( 2QQ2 ). (Radiative cooling of the thermal material might be ex- 
pected to enhance the density contrast here still further.) Structures 
like this, positioned roughly inside the boundary of the radio lobes 
(though note that projection will affect the lobes and the bar/arms 
structure differently, so that we do not expect an exact match be- 
tween the two) are a smoking gun for a system in which the lobes 
have come into pressure balance at large radii and are now being 
squeezed out as in Model C of |Scheuer| ( | 1 974| ) . In real observations, 
structures that could be interpreted in this way are in fact ubiquitous 
in observations of FRIIs; examples include the central, bright, cool 
structure in Cygnus A (jWilson et al.|r2006 ), the very obvious bar 
and arms in 3 C444 (|Croston et al.||2011| ), the 'ridges' in 3C285 
and 3C442A pardcastle et aL[ |2007| ) and in 3C401 ( |Reynolds| 
|et al.| |2005| ). In the best-studied cases with the clearest evidence 
for ongoing strong shocks (Cygnus A and 3C444) the 'bar' and 
'arm' temperatures are consistent with the unshocked cluster tem- 
perature, precisely as expected in our models. The agreement be- 
tween our simulations and the observations of 3C 444, which is an 
FRII at z = 0.153 with a slightly higher jet power embedded in 
a slightly richer environment, is very striking; forthcoming deep 
XMM-Newton observations of this object may provide the opportu- 
nity for a detailed test of the model's temperature predictions. This 
should be combined with sensitive low-frequency observations, e.g. 
with LOFAR, to establish exactly how the lobes relate to the X-ray 
structures. 



4.6 Suppression of cooling? 

Our simulation setup was explicitly designed so that the jet power 
of the FRII exceeds by roughly an order of magnitude the bolomet- 
ric X-ray luminosity of the simulated cluster environment. But this 
is not a sufficient (nor indeed a necessary) condition for suppressing 
radiative cooling of the X-ray-emitting material, because the cool- 
ing rate is a strong function of density and thus of radius, and the 
key requirement for prevention of the traditional cooling catastro- 
phe is to suppress cooling at the centre of the group/cluster. While 
numerical modelling shows that powerful, FRII-like radio galaxies 

l[2003j), it has been 



can suppress cooling (e. g.|Basson & Alexander] 
clear for some time (e.g. jOmma & Binney 2004| ) that they are not 



an efficient way of doing so, because the bulk of the energy /entropy 
injection takes place on scales that are much larger than the cool- 
ing radius (i.e. the radius within which the cooling time is less than 
some threshold value). In sources like ours, which are intended to 
come into pressure balance at the centre at late times, the situation 
is even worse, because there is little or no energy input at the centre 
of the cluster at all, while the increase in axial ratio with length also 
means that less of the available volume of the cluster is affected at 



late times. However, the cluster cores do undergo significant heat- 
ing early in the simulations. 

We can quantify these effects in our simulations by looking 
at radial profiles of some quantities related to cooling as a func- 
tion of simulation time. Fig. [TT] shows an example provided by the 
M25-35-40 run, but all our results are broadly similar. What we 
see is that, as expected, the total energy input to the atmosphere is 
mostly co-located with the front of the lobes. The energy input to 
the atmosphere at small radii is actually negative, because the ef- 
fect of the radio galaxy is to remove gas from the central regions 
without providing much in the way of long-term heating in these 
locations. This can be seen in the density profiles in Fig.[T7]- in the 
short term, the radio source has a dramatic impact on the density 
of hot gas in the central regions, but as the central pressure drops 
and material flows in behind the radio lobes the mean density pro- 
file rapidly recovers to something not very dis- similar to its original 
shape behind the shock front. Similarly, the entr opy profiles - here 



we are plotting the entropy index, a — Pn |Kaiser & Binney 



|2003| ) in simulation units - show a very marked effect at early times 
but only a very modest increase in the entropy in the central regions 
of the plot at late times. The effect on cooling time (proportional to 
T^/^/n: not plotted) is very similar to that on entropy. Note also 
the almost completely negligible effect of the radio source on the 
entropy profile at large radii - this is because these are volume- 
averaged radial profiles, and at these large radii very little of the 
volume of the cluster is actually affected by the radio galaxy. 

Why do we see this behaviour of the entropy index? In our 
simulations, entropy is generated at shock fronts, redistributed in 
mixing regions, and conserved otherwise. The strong entropy in- 
crease in the whole core region at early times means that this gas 
must in the end rise - either buoyantly or by being carried out by the 
expanding radio source - out of the centre of the cluster; in other 
words, much of the core has been efficiently heated. However, the 
core is then refilled by adjacent gas from larger initial radii. Cor- 
respondingly, the density decreases in this region. Therefore, while 
the effect of the jet was to increase the entropy of much of the gas 
that was in the core initially by a factor of several tens, the gas that 
is in the core at the end of the simulation has only a 10-20 per cent 
higher entropy than the initial conditions in this region. Because the 
radio galaxies continue to drive strong shocks at large radii, they 
continue to increase the entropy of the shocked gas at all times: this 
can be seen in a plot of the integrated entropy, / adN, as a function 
of lobe length (Fig.[T8]). However, as the profiles show, this entropy 
increase is no longer relevant to the cooling core, or significant over 
the total volume of the cluster. 

We conclude that long-lived FRII radio sources with realistic 
dynamics are not expected to be very efficient at suppressing cool- 
ing in this type of environment. Though much of the energy is ther- 
malised and ends up in the ambient gas, only the central parts of the 
cluster are efficiently, irreversibly heated, and much of this material 
then leaves the central regions; at late times (i.e. '-^ 10^ years after 
the radio activity begins), there are only very modest effects on the 
central entropy and cooling time. To couple AGN output effectively 
to the gas with short cooling times in the centres of clusters, radio 
sources need not to grow beyond the central cluster scale, so that 
they are continually in the early-time regime seen in Fig. [T7] An 
excellent way of doing this is to make the radio source intermittent 
(see e.g. |Mendygral et aL||2012| and observations of the centres of 
cool-core clusters such as Perseus show that this is what happens in 
nature most - but not all! - of the time. One of the key questions to 
be solved in 'AGN feedback' models must be how the AGN in cool 
cores 'know' that they are required to be intermittent in this way 
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Figure 17. Radial profiles of excess energy, number density and entropy in the external medium as a function of radius for 6 different times in the M25 -35-40 
simulation. The excess energy profile (note linear scale on the y axis) shows total energy (thermal and kinetic) in all the cells that are not inside the lobe at a 
given radius, minus the value in those cells at t = 0. The density and entropy profiles are the volume-weighted mean values at each radius. See Figs[^and[2| 
for images of the simulated density and temperature in this simulation at these times (t > 0). 
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Figure 18. Integrated excess entropy as a function of lobe length for the simulations with = 25. A baseline corresponding to the initial (t = 0) entropy of 
the gas inside the shock front has been subtracted from all values plotted. The slope of this plot gives a measure of the efficiency with which the radio galaxy 
irreversibly heats the external medium as it grows. The three distinct regions here can be understood in terms of the three phases of lobe growth speed, and 
therefore shock speed (Fig.[3j; weaker shocks have a lesser effect on entropy ^jKaiser & Binney|f2003) . 



while other radio galaxies, like the classical doubles, can exhibit 
long-term persistent jet activity. 



5 SUMMARY AND CONCLUSIONS 

We have carried out high-resolution, 2D axisymmetric hydrody- 
namical simulations of the evolution of the lobes of radio galaxies 
in a wide range of realistic environments. The simulations are in- 
tended to be closely matched to the physical conditions known to 
obtain in the environments of FRII radio galaxies, and to probe all 
size scales from the collimation of an initially conical jet up to (and 
beyond) the size scale at which the lobe pressure is expected to drop 
below the pressure in the external medium. Although the details of 
the dynamics internal to the lobe are probably not realistic, we ar- 
gue that the energetics of the interaction with the external medium 
should be unaffected by this so long as there is a mechanism that 
efficiently 'thermalizes' the kinetic energy supplied by the jet, as 
is the case in our simulations once the lobes grow to a few tens of 
simulation lengths in size. 

The key results of our work appear to be independent of both 
environment and initial jet Mach number, within the range that we 
have been able to probe, and can be summarized as follows: 

(i) As expected both from observations and from the design of 
our simulations, the lobes do come into pressure balance with the 
external medium in the course of the simulations (Section [4T| Sec- 
tion |4.2| ). At late times, the lobes are driven away from the centre of 
the host cluster by dense, high-pressure ICM. Thus the lobe growth 
is more consistent with Model C than with Model A of Scheuer 
(1974), or, e quivalently, we observe the 'cocoon crushing' process 
of Williams] { [ 1 99 1 1 . Lobe growth is not self-similar, and indeed ax- 



ial ratios tend to increase with time on hundred-kpc scales, consis- 
tent with what is observed (Section [43] ). 

(ii) Unsurprisingly, given the above, the models of KA do not 
describe some of the behaviour of the simulated radio galaxies 
(Section [43] l. The late-time length of the sources as a function of 
time is not badly described by the relations derived by KA, which 
is not surprising since the environment has close to a power-law 
behaviour of density with distance on these scales and the momen- 
tum flux of the jet, distributed over the 'working surface', continues 
to drive the linear expansion of the radio source. However, the KA 
relations significantly overestimate the growth of lobe volume with 
time. 

(iii) The pressure of the undisturbed external environment at the 
midpoint of the lobes is a reasonable, though not perfect, estimator 
of the pressure in the lobes once normal lobe growth is established 
(Section [43] ). Lobes remain strongly overpressured with respect to 
the ambient medium at their tips. This reproduces a key observa- 
tional result based on estimates of the lobe pressure from inverse- 
Compton observations. 

(iv) The work done by the radio galaxy on the environment (in- 
cluding the thermal, kinetic and potential energy of the shocked 
gas) is very close to the internal energy (thermal and kinetic) of 
the lobes at all times once the lobes are established (Section [43] ). 
A simple pextV estimate of the work done by the lobes underesti- 
mates it by a factor of a few (up to an order of magnitude in extreme 
cases). 

(v) The simulated radio galaxies can reproduce the approximate 
radio luminosities of real 3C radio sources whose jet powers are 
believed to be comparable (Section [4!4| . However, there is a large 
scatter (well over an order of magnitude) in the radio powers ex- 
pected from sources of the same power (as all our simulated sources 
are) in different environments. 
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(vi) X-ray visualization of our simulations (Section [43] ) makes 
some predictions for the morphological and temperature structure 
that we should observe that are in principle testable, and that show 
striking similarities to what is observed in a few well-studied FRII 
sources. 

(vii) Finally, we note that, as has been known for a while, FRII 
radio galaxies with realistic dynamics do not couple effectively to 
the hot medium at the centres of clusters, and therefore are not 
in general effective at suppressing cooling (Section though 
their energetic impact can be very considerable in localized regions 
around the lobes in the outer parts of the cluster. 
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Several of these points have importance for observational 
work and for models of the energetic impact of radio galaxies on 
their environments. Firstly, our conclusion that the internal and ex- 
ternal energy of the radio galaxy remain in rough balance, indepen- 
dent of the properties of the external environment or the initial jet 
Mach number, means that there is a very easy recipe for estimating 
the work done by an FRII radio galaxy on its environment; simply 
estimate the energy stored in the lobes, ideally using an inverse- 
Compton observation to measure the lobe magnetic field strength, 
and one also has, to within better than a factor 2, the work done 
on the external environment. Although it might be objected that 
our modelling does not include important features of the physics of 
radio galaxies, such as magnetic fields, which might be expected 
to change this conclusion, we are reassured by the similar conclu- 
sions reached for specific simulations with more realistic physics, 
including magnetic fields, by jO'Neill et al.| ( [2005] ) and |Gaibler et al.| 
( [20091 ). 

Secondly, if we accept that all FRII radio galaxies follow the 
basic pattern seen in our modelling and come into rough pressure 
balance with their environments, independent of the initial power 
of their jets, at around the lobe midpoint, once they have grown 
to a sufficiently large scale (i.e. once they are in a roughly power- 
law atmosphere?), then this means that large-scale FRIIs provide a 
probe of the pressure of their environments; if one observes an FRII 
on hundred-kpc scales, the external pressure at its midpoint should 
be of order the internal pressure, which can again be assessed by 
inverse-Compton observations or, with some assumptions about 
magnetic field strength, from synchrotron emission alone. In the 
coming decades, where deep radio observations may be far easier 
to get than deep X-ray observations, FRII radio lobes may become 
a useful proxy for high-redshift, high-pressure environments. 

Thirdly, and on a less positive note, we draw attention to 
the very large scatter seen in the radio luminosity plots. Although 
we have not included radiative losses (synchrotron and inverse- 
Compton) in our modelling, these are unlikely to decrease the scat- 
ter, and indeed observational constraints on the jet power/radio 
power relationship for FRIIs show similar amounts of scatter (e.g. 
[Daly et aL] [2012 ). While there should be a statistical correlation 
between jet power and radio power, our results suggest that naive 
application of the results of, e.g., |Willott et al.| ( |l999| ), while they 
may be adequate in a statistical sense, may give very significant 
errors for individual radio sources. 

In future work, we hope to check the robustness of the conclu- 
sions we have drawn here by including more of the relevant physics 
in the simulations (e.g. more detailed 3D simulations, MHD, radia- 
tive losses...) while continuing to match our modelling explicitly to 
realistic radio galaxies and their environments. 
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